Years

year_from <- min(bfs_web_deaths_monthly$Year)

year_reg <- year_from + 10

year_max <- 2020

Pandemics

flu_years_hc <- c(1890, 1918, 2020)

flu_years_hc_affected <- c(seq(1890 + 1, 1890 + 10),
                           seq(1918 + 1, 1918 + 10))

flu_years <- c(1890, 1918, 1957, 1968, 1977, 2009, 2020)

flu_years_affected <- c(flu_years_hc_affected,
                        seq(1957 + 1, 1957 + 10),
                        seq(1968 + 1, 1968 + 10),
                        seq(1977 + 1, 1977 + 10),
                        seq(2009 + 1, 2009 + 10))

Pops

Data till 2019 - taking mean from jan & dec values.

bfs_web_pop_yearly <- read_rds("data/BfS/bfs_web_pop_yearly.Rds") %>% 
  filter(Year >= min(bfs_web_deaths_monthly$Year)) %>% 
  mutate(pop = round((pop_jan + pop_dec) / 2)) %>% 
  select(-pop_jan, -pop_dec)

2020 simply predicted using pops from last 10 years

(pop_2020_pred <- predict(
  lm(pop ~ Year, 
     data = bfs_web_pop_yearly %>% filter(Year >= 2010)),
  new = tibble(Year = 2020)))

      1 
8693850 

Together

bfs_web_deaths_monthly %<>% 
  left_join(bfs_web_pop_yearly) %>%
  mutate(flu_year_hc = ifelse(Year %in% flu_years_hc, 1, 0),
         flu_year_hc_affected = ifelse(Year %in% flu_years_hc_affected, 1, 0),
         deaths_flu_year_hc = ifelse(Year %in% flu_years_hc, NA, deaths),
         flu_year = ifelse(Year %in% flu_years, 1, 0),
         flu_year_affected = ifelse(Year %in% flu_years_affected, 1, 0),
         deaths_flu_year = ifelse(Year %in% flu_years, NA, deaths),
         pop = ifelse(Year == 2020, pop_2020_pred, pop),
         si_one = sin(2*pi*Month/12),
         si_two = sin(4*pi*Month/12),
         co_one = cos(2*pi*Month/12),
         co_two = cos(4*pi*Month/12)) 

rm(bfs_web_pop_yearly)

4 harmonic variables

bfs_web_deaths_monthly %>% 
  filter(Year >= 2018) %>% 
  select(si_one, co_one) %>% 
  mutate(row = row_number()) %>% 
  pivot_longer(-row) %>% 
  ggplot() +
  geom_line(aes(row, value, color = name)) +
  xlab("") +
  theme_minimal()

bfs_web_deaths_monthly %>% 
  filter(Year >= 2018) %>% 
  select(si_two, co_two) %>% 
  mutate(row = row_number()) %>% 
  pivot_longer(-row) %>% 
  ggplot() +
  geom_line(aes(row, value, color = name)) +
  xlab("") +
  theme_minimal()

Three pandemics

With denominator

5- & 10-year medians

bfs_web_deaths_monthly %<>% 
  arrange(Month, Year) %>% 
  mutate(median =             slide_dbl(deaths, 
                                        stats::median, na.rm = TRUE, 
                                        .before = 10, .after = -1, .complete = TRUE),
         median_flu_year_hc = slide_dbl(deaths_flu_year_hc, 
                                        stats::median, na.rm = TRUE,
                                        .before = 10, .after = -1, .complete = TRUE),
         median_flu_year =    slide_dbl(deaths_flu_year, 
                                        stats::median, na.rm = TRUE, 
                                        .before = 10, .after = -1, .complete = TRUE),
  ) %>% 
  arrange(Year, Month)

Impact of 1957 on 1958

Seasonality

dcmp <- bfs_web_deaths_monthly %>%
  mutate(Time = row_number()) %>% 
  tsibble::as_tsibble(index = Time) %>% 
  fabletools::model(feasts::STL((deaths/pop)*100000 ~ season(window = Inf)))

# fabletools::components(dcmp)

fabletools::components(dcmp) %>% 
  autoplot()

rm(dcmp)

Outcome

CI via bootstrap

Method from https://statcompute.wordpress.com/2015/12/20/prediction-intervals-for-poisson-regression/

p_load(doParallel, foreach)

registerDoParallel(cores = 4)

boot_pi <- function(model, pdata, n, p) {
  odata <- model$data
  lp <- (1 - p) / 2
  up <- 1 - lp
  seeds <- round(runif(n, 1, 1000), 0)
  boot_y <- foreach(i = 1:n, .combine = rbind) %dopar% {
    set.seed(seeds[i])
    bdata <- odata[sample(seq(nrow(odata)), size = nrow(odata), replace = TRUE), ]
    bpred <- predict(update(model, data = bdata), type = "response", newdata = pdata)
    rpois(length(bpred), lambda = bpred)
  }
  boot_ci <- t(apply(boot_y, 2, quantile, c(lp, up)))
  return(data.frame(pred = predict(model, newdata = pdata, type = "response"), 
                    lower = boot_ci[, 1], upper = boot_ci[, 2]))
}

Looping over years

Analyses using 10 previous years.
With three different inclusions depending on pandemic / flu years

for (YEAR in year_reg:year_max){
  
  print(paste("Analysing year", YEAR))
  
  reg_data <- bfs_web_deaths_monthly %>% 
    filter(Year >= YEAR - 10 & Year < YEAR) 
  
  pred_data <- bfs_web_deaths_monthly %>% 
    filter(Year == YEAR) %>%
    mutate(
      fit = NA_real_,
      lpi = NA_real_,
      upi = NA_real_,
      fit_flu_hc = NA_real_,
      lpi_flu_hc = NA_real_,
      upi_flu_hc = NA_real_,
      fit_flu = NA_real_,
      lpi_flu = NA_real_,
      upi_flu = NA_real_)
  
  # summary(reg_data$Year)
  # summary(pred_data$Year)
  
  # all data
  summary(monthly <- glm(deaths ~ Year + 
                           si_one + si_two + co_one + co_two, 
                         offset = log(pop),
                         data = reg_data, 
                         family = "poisson"))
  
  predict <- boot_pi(monthly, pred_data, 1000, 0.99)
  
  pred_data %<>%
    mutate(
      fit = predict$pred,
      lpi = predict$lower,
      upi = predict$upper,
      fit_flu_hc = predict$pred,
      lpi_flu_hc = predict$lower,
      upi_flu_hc = predict$upper,
      fit_flu = predict$pred,
      lpi_flu = predict$lower,
      upi_flu = predict$upper
    )
  
  # flu hc excluded
  if (YEAR %in% flu_years_hc_affected) {
    
    print(paste("    >> Extra analyses for hc flu year exclusions"))
    
    summary(monthly_flu_hc <- glm(deaths_flu_year_hc ~ Year + 
                                    si_one + si_two + co_one + co_two, 
                                  offset = log(pop),
                                  data = reg_data, 
                                  family = "poisson"))
    
    predict_flu_hc <- boot_pi(monthly_flu_hc, pred_data, 1000, 0.99)
    
    pred_data %<>%
      mutate(
        fit_flu_hc = predict_flu_hc$pred,
        lpi_flu_hc = predict_flu_hc$lower,
        upi_flu_hc = predict_flu_hc$upper
      )
  }
  
  # all flu excluded
  if (YEAR %in% flu_years_affected) {
    
    print(paste("    >> Extra analyses for flu year exclusions"))
    
    summary(monthly_flu <- glm(deaths_flu_year ~ Year + 
                                 si_one + si_two + co_one + co_two, 
                               offset = log(pop),
                               data = reg_data, 
                               family = "poisson"))
    
    predict_flu <- boot_pi(monthly_flu, pred_data, 1000, 0.99)
    
    pred_data %<>%
      mutate(
        fit_flu = predict_flu$pred,
        lpi_flu = predict_flu$lower,
        upi_flu = predict_flu$upper
      )
  } 
  
  if (YEAR == year_reg) {
    expected_deaths_monthly <- pred_data
  } else {
    expected_deaths_monthly <- bind_rows(expected_deaths_monthly, pred_data)
  }
}

[1] "Analysing year 1887"
[1] "Analysing year 1888"
[1] "Analysing year 1889"
[1] "Analysing year 1890"
[1] "Analysing year 1891"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1892"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1893"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1894"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1895"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1896"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1897"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1898"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1899"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1900"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1901"
[1] "Analysing year 1902"
[1] "Analysing year 1903"
[1] "Analysing year 1904"
[1] "Analysing year 1905"
[1] "Analysing year 1906"
[1] "Analysing year 1907"
[1] "Analysing year 1908"
[1] "Analysing year 1909"
[1] "Analysing year 1910"
[1] "Analysing year 1911"
[1] "Analysing year 1912"
[1] "Analysing year 1913"
[1] "Analysing year 1914"
[1] "Analysing year 1915"
[1] "Analysing year 1916"
[1] "Analysing year 1917"
[1] "Analysing year 1918"
[1] "Analysing year 1919"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1920"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1921"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1922"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1923"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1924"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1925"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1926"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1927"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1928"
[1] "    >> Extra analyses for hc flu year exclusions"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1929"
[1] "Analysing year 1930"
[1] "Analysing year 1931"
[1] "Analysing year 1932"
[1] "Analysing year 1933"
[1] "Analysing year 1934"
[1] "Analysing year 1935"
[1] "Analysing year 1936"
[1] "Analysing year 1937"
[1] "Analysing year 1938"
[1] "Analysing year 1939"
[1] "Analysing year 1940"
[1] "Analysing year 1941"
[1] "Analysing year 1942"
[1] "Analysing year 1943"
[1] "Analysing year 1944"
[1] "Analysing year 1945"
[1] "Analysing year 1946"
[1] "Analysing year 1947"
[1] "Analysing year 1948"
[1] "Analysing year 1949"
[1] "Analysing year 1950"
[1] "Analysing year 1951"
[1] "Analysing year 1952"
[1] "Analysing year 1953"
[1] "Analysing year 1954"
[1] "Analysing year 1955"
[1] "Analysing year 1956"
[1] "Analysing year 1957"
[1] "Analysing year 1958"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1959"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1960"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1961"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1962"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1963"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1964"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1965"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1966"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1967"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1968"
[1] "Analysing year 1969"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1970"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1971"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1972"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1973"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1974"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1975"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1976"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1977"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1978"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1979"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1980"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1981"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1982"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1983"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1984"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1985"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1986"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1987"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 1988"
[1] "Analysing year 1989"
[1] "Analysing year 1990"
[1] "Analysing year 1991"
[1] "Analysing year 1992"
[1] "Analysing year 1993"
[1] "Analysing year 1994"
[1] "Analysing year 1995"
[1] "Analysing year 1996"
[1] "Analysing year 1997"
[1] "Analysing year 1998"
[1] "Analysing year 1999"
[1] "Analysing year 2000"
[1] "Analysing year 2001"
[1] "Analysing year 2002"
[1] "Analysing year 2003"
[1] "Analysing year 2004"
[1] "Analysing year 2005"
[1] "Analysing year 2006"
[1] "Analysing year 2007"
[1] "Analysing year 2008"
[1] "Analysing year 2009"
[1] "Analysing year 2010"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2011"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2012"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2013"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2014"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2015"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2016"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2017"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2018"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2019"
[1] "    >> Extra analyses for flu year exclusions"
[1] "Analysing year 2020"

rm(predict, predict_flu_hc, predict_flu, pred_data, reg_data, YEAR)
write_rds(expected_deaths_monthly, "data/expected_deaths_monthly.Rds")

Two predictions

Animate

Stacked

Diffs